Εργαστηριακή Άσκηση 2#
Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.
Live Code
Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).
Μέρος 1: Εισαγωγή#
Στο MATLAB οι συναρτήσεις \(fft\) και \(ifft\) υποθέτουν ζεύγος μετασχηματισμού Fourier \(x(t)\) και \(X(f)\) υπολογισμένων σε μη αρνητικά διαστήματα \( t=[0:N-1]ts \) και \( f=[0:N-1]fo. \) Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το \(X(f)\) με τη βοήθεια της συνάρτησης \(fft\) (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα \(x(t)\) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της \(ifft\).
Το MATLAB διαθέτει τη συνάρτηση \(fftshift\) για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) ή στις συχνότητες \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). Με τον τρόπο αυτό μπορούμε να παράγουμε τα \(xb(t)\) και \(Xb(f)\) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier.
Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα \([1 2 3 4]\) ως το αποτέλεσμα του \(FFT\) μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας \(fs/2\), που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην \(–fs/2\) είτε στην \(+fs/2\). Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες \(+fs/4\) και \(–fs/4\). Εφαρμόζοντας την \(fftshift\), το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα \(–fs/2\), το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα \(–fs/4\) ακολουθούμενο από το dc και τη συχνότητα \(+fs/4\). Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το \(±fs/2\). Έτσι για το διάνυσμα \([1 2 3]\), η εφαρμογή της \(fftshift\) θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες \(–fs/3, 0, +fs/3\).
Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις \(fft\) και \(ifft\) αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες
\(h == \text{ifft}(\text{fft}(h))\) και \(H == \text{fft}(\text{ifft}(H))\)
Η πρώτη υποδεικνύει ότι η είσοδος της \(ifft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(fft\), και η δεύτερη ότι η είσοδος της \(fft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(ifft\).
Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import Layout
from IPython.display import display, clear_output
np.random.seed(0)
# Generate random data
t = np.linspace(-40, 40, 80) # 80 samples
random_data = np.random.randn(len(t))
# Initial start and end points for the window
start, end = 30, 50
# Output widget for the plots
output = widgets.Output()
def apply_window_and_plot(start, end):
with output:
clear_output(wait=True)
# Create a sinusoidal window
sin_window = np.zeros_like(t)
sin_window[start:end] = np.sin(np.linspace(0, np.pi, end - start))
# Apply sinusoidal window to random data
windowed_data = sin_window * random_data
# Append the negative times to the end of the time array
t_append_negative = np.concatenate((t[t >= 0], t[t < 0] + t[-1] - t[0] + (t[1] - t[0])))
windowed_data_append_negative = np.concatenate((windowed_data[t >= 0], windowed_data[t < 0]))
plt.close('all') # Close all existing figures
# Create subplots
fig, axs = plt.subplots(2, 1, figsize=(10, 10))
# Plot original windowed data
axs[0].stem(t, windowed_data, linefmt='C0-', markerfmt='C0o', basefmt=" ")
axs[0].set_title("Original Sinusoidal Windowed Random Data", fontsize=16)
axs[0].set_xlabel("Time (samples)", fontsize=12)
axs[0].set_ylabel("Amplitude", fontsize=12)
# Plot with negative times appended
axs[1].stem(t_append_negative, windowed_data_append_negative, linefmt='C1-', markerfmt='C1o', basefmt=" ")
axs[1].set_title("Sinusoidal Windowed Random Data with Negative Times Appended", fontsize=16)
axs[1].set_xlabel("Time (samples)", fontsize=12)
axs[1].set_ylabel("Amplitude", fontsize=12)
# Adjust layout and show plot
plt.tight_layout()
plt.show()
# Define the range slider
window_range_slider = widgets.IntRangeSlider(
value=[start, end],
min=0,
max=len(t),
step=1,
description='Window Range:',
layout=Layout(width='90%'), # Setting the width to 90%
continuous_update=False
)
# Function to update the plot based on the range slider
def on_range_change(change):
new_start, new_end = change['new']
apply_window_and_plot(new_start, new_end)
# Observe the range slider for changes
window_range_slider.observe(on_range_change, names='value')
# Group the slider and output widget in a vertical box layout
vbox = widgets.VBox([window_range_slider, output])
# Display the VBox
display(vbox)
# Initialize with the default plot
apply_window_and_plot(start, end)
